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ABSTRACT 

Within the standard models of particle physics and cosmology we have calculated the big-bang 
prediction for the primordial abundance of ^He to a theoretical uncertainty of less than 0.1% 
{5Yp < ±0.0002), improving the current theoretical precision by a factor of 10. At this accuracy 
the uncertainty in the abundance is dominated by the experimental uncertainty in the neutron mean 
lifetime, r„ = 885.4ib2.0sec. The following physical effects were included in the calculation: the zero 
and finite-temperature radiative, Coulomb and finite- nucleon-mass corrections to the weak rates; 
order-a quantum-electrodynamic correction to the plasma density, electron mass, and neutrino 
temperature; and incomplete neutrino decoupling. New results for the finite-temperature radiative 
correction and the QED plasma correction were used. In addition, we wrote a new and independent 
nucleosynthesis code designed to control numerical errors to be less than 0.1%. Our predictions 
for the ^He abundance are presented in the form of an accurate fitting formula. Summarizing 
our work in one number, Yp{r] = 5 x 10"^°) = 0.2462 ± 0.0004 (expt) ± < 0.0002 (theory). 
Further, the baryon density inferred from the Buries- Tytler determination of the primordial D 
abundance, Qsh^ = 0.019 ± 0.001, leads to the prediction: Yp = 0.2464 ± 0.0005 (D/H) ± < 
0.0002 (theory) it 0.0005 (expt) . This "prediction" and an accurate measurement of the primeval 
^He abundance will allow an important consistency test of primordial nucleosynthesis. 



1 Introduction 



Big-bang nucleosynthesis (BBN) is one of the observational pillars of the standard cosmology. 
Further, it has the potential to be a precision probe of the early universe and fundamental physics Q, 
Observations of light-element abundances have improved dramatically over the past few years, 
and the current and planned precision measurements of D, '^He, ^He and ^Li, should allow a precise 
(10% or better) determination of the baryon density and consistency check of BBN, but only if 
the theoretical predictions of the light element abundances are as good as the observations. In 
particular, a measurement of the primeval D abundance pins down the baryon density, and in turn 
makes predictions for the other three abundances. Because the subsequent evolution of the ^He 
abundance is simple - stars make ^He- and measurements have the potential of determining Yp to 
three significant figures 0, ^, ^, 0, ^, 10|, ^He can provide an important consistency check of 
BBN. Furthermore, an independent determination of the baryon density from cosmic microwave 
background anisotropics will soon test the consistency of the standard model of cosmology. Finally, 
the combination of accurate observations and theory can be used to test physics beyond the standard 
model of particle physics e.g., by imposing a strict limit on the number of light neutrino 

species Cosmology is entering a high precision age, and this motivates high-precision 

BBN predictions. 

Over the years, theoretical study of "^He synthesis has been intense, with the following effects 
being considered: Coulomb and radiative corrections to the weak rates |15, 16, 17, 18, 19, |2^, BBN 
code numerical errors |18], nuclear reaction rate uncertainties |21, 22 1, finite-temperature QED 
plasma corrections [15, 23 1, the effect of finite-nucleon mass [24, 25, 26 1, and incomplete neutrino 
decoupling [15, However, the corrections have been incorporated in a patchwork fashion and 



a recent informal poll of BBN codes indicated a spread of 1 % in the predicted value of the ^He 
abundance for fixed ry and r„. 

The goal of this work was a calculation of the primordial abundance of ^He, within the standard 
models of particle physics and cosmology, accurate enough so that its uncertainty is dominated by 
the experimental uncertainty in the neutron mean lifetime^, = 885.4 ± 2.0 sec [M, E9L ROi. 
Because r„ is so accurately known ((5r„/T„ = 0.23%), it is used to normalize all of the weak rates 
that interconvert neutrons and protons: ep <-> z^n, e^n <-> and n <-> pev'. The baryon-number 
fraction of ^He produced (= Yp) depends sensitively on the weak rates because they determine the 
neutron-to-proton ratio n/p before nucleosynthesis, and essentially all of the neutrons around at the 
onset of nucleosynthesis go into forming '^He. We have determined the effect on Yp by perturbing 
the weak rates in the standard code [31|, 



6Yp 



(1) 



Since the weak rates scale as l/x„, this estimate implies that introduces an uncertainty in Yp 
of 0.18%. We use this uncertainty to set our goal for all theoretical uncertainty. 



^The Particle Data Group currently recommends r„ = 887 ± 2 sec [^g]- ^ recent measurement using ultracold 
neutrons indicates a slightly lower value, r„ — 885.4 ±0.9 (stat) ±0.4 (sys) sec |^^. For our central value we use 885.4 
sec and for the uncertainty we use ±2 sec. 
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To meet our goal we need to calculate the weak rates to precision of better than 0.23%. Another 
source of errors in Yp come from thermodynamics, i.e., the energy density p, the pressure P and the 
neutrino temperature T^. To determine how accurately we need to know thermodynamic quantities, 
we can estimate the change in Yp due to a change in a thermodynamic quantity, e.g., p. Again, 
using the standard code, we find 

f£==0.4^, (2) 

Yp p 

This indicates that we should calculate thermodynamic quantities to better than 0.45%. 
When calculating Yp to this precision, several factors must be considered: 

1. Weak rate and thermodynamics numerics: most quantities to be calculated involve integra- 
tions that must be done numerically. 

2. ODE integration numerics: nucleosynthesis codes contain finite stepsize errors. 

3. Nuclear reaction rates: errors originate from experimental uncertainties in the nuclear reaction 
data, as well as from neglecting nuclear reactions important to BBN. 

4. Weak-rate physics: there are several small physical effects that must be calculated, including 
Coulomb, zero and finite-temperature radiative corrections, and the effect of finite- nucleon 
mass. 

5. Thermodynamics physics: for temperatures much greater than the electron mass, there are 
order-a quantum electro dynamic corrections to the equation of state of the plasma. 

6. Incomplete neutrino decoupling: neutrinos share partially in the entropy release when 
pairs annihilate. 

Items I, I and | are addressed in the next Section; item ^ is addressed in Sec. 3. Items |5| and ^ 
are taken up in Sec. 4, and a summary of our results is given in the final section. 

We mention that we have not considered the 0(a^/^) collective plasma effects due to the presence 
of the copious numbers of pairs at the time of BBN, because they are safely below our theoretical 



error budget of 0.1 % for Yp. These effects, all of relative size 0.1 % and calculated in Ref. |Q, are: 
the enhancement of nuclear reaction rates due to Debye screening of nuclear charge; the contribution 
of longitudinal plasmon modes (A: < ~ 47rng±/r) to the energy density and pressure; the 
(negative) contribution to the energy density and pressure of the electromagnetic interaction of 
pairs; and the reduction of the energy and pressure of photons due to plasma effects on low 
frequency photons {k < LOp). Finally, while we have tried to be exhaustive and very careful in 
our analysis, we cannot rule out systematic theoretical error: that is, the possibility that we have 
neglected some microphysical effect as important as those we have included. 

2 Numerics 
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2.1 BBN Code 



We have written a new nucleosynthesis code that is independent of the standard (Kawano) code |31|. 
The heart of any nucleosynthesis code is the set of ordinary differential equations that govern the 
evolution of the abundances of the light elements (see, e.g., Ref. Q). Our code tracks protons, 
neutrons, D, T, ^He, ^He, ^Li, ''Li and ''Be. The baryon-number fraction of element i is given by|^ 



ns 1 + I]i {ni/nn) ' 

where Ai is the element's atomic number, ni it's number density, and is the baryon-number 
density. (Note, by convention Yp is used to denote X4). Nuclear reaction rates govern the evolution 
of the elemental abundances. Conservation of baryon number provides the constraint: 



1.000. (4) 



We take for our initial temperature, Tj = 10 MeV, and for our initial abundances, the nuclear 
statistical equilibrium (NSE) values: 

Xa = 9A rc(3)^"ivr(^-^)/22(3^-5)/2l ^5/2 ( ^\ r/^-ixf X;^-^e^^/^ (5) 

L J yniN J ^ 

where A is the atomic number, niN — 940 MeV is the nuclear mass, rj is the baryon-to-photon ratio, 
Ba is the binding energy of species A, and C(3) ^ 1.20206. At temperatures greater than about 1 
MeV, the nuclear rates are sufficiently high to cause the abundances to rapidly assume their NSE 



values. ( As discussed in Ref |35|, the final abundances are very insensitive to the assumed initial 

abundances). If we make the well-justified assumption that the elements are always in kinetic 

equilibrium, then the rate coefficients depend only on r] and T. This implies the important and 

well known conclusion that the predictions of nucleosynthesis are a function of only one parameter, 

rj, which is equivalent to since = 2.7277 it 0.002 K is so well known. 

Several important quantities enter into the evolution equations: weak rates, thermodynamic 

quantities and nuclear reaction rates. For the weak rates, we define the total conversion rates (per 

^Baryon-number fraction and baryon-mass fraction differ by order 1% due to nuclear binding energy. Because 
nuclear reactions change the total mass in baryons, the mass fraction of species Ai (= Xl^"""") can change even if the 
number of species Ai does not. The mass fraction of species Ai is 

nirrii m rrii 1 



XT 



where rrii is the mass of species i: e.g., 7714 = 4.002602 amu and niH ~ 1.00783 amu. For Yp — 0.25 and the primordial 
mix of elements X™""" = 0.24866. Similarly, the relationship between the baryon mass density and rj depends on 
elemental composition. For the primordial mix with Yp — 0.25, 

Qb = 3.66 X 10'' 7?, 

with Tj = 2.7277 K. Assuming a mass of 1 amu per nucleon, the prefactor is 3.639 x 10*^, and for solar abundance, 
the prefactor is 3.66043 x 10^ 
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Table 1: Reactions used in our code. 

neutron or proton): 

Tp^n = ^ep—*i/n ~l~ ^up—yc^n ~^ ^pev^n- (6) 

Simple expressions for these rates may be obtained assuming no radiative corrections and infinite 
nucleon mass. The thermodynamic quantities that must be calculated are p{T), T^(T), pb{T) and 
the differential time-temperature relation dt/dT. 

Our BBN code is independent from the standard code, with one exception: It uses the same 
nuclear-rate data (with the exception of the weak rates). The nuclear-reaction network corresponds 
to the smallest one offered by the standard code, which contains the reactions listed in Table |l|. 
Although this network is much smaller than the largest offered in the standard code, we have verified 
that the effect on Yp of neglecting these additional reactions is less than 10~^. The light-element 
abundances predicted by our code are shown in Fig. |^. 

2.2 Numerical Accuracy of the BBN Codes 

Because the differential equations governing the light-element abundances are stiff, an implicit 
integrator was used to evolve them. Instead of specifying explicit time steps, as in the standard code, 
the desired final accuracies are specified as parameters of our code's integrator. The temperature 
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steps are then determined adaptively. Integrator accuracy parameters are chosen to be smah enough 
so that stepsize errors were much smaher than the allowed error in Yp. 

To calculate the weak rates and thermodynamic quantities accurately, we proceed as follows 
(see, e.g., Ref. ||3^]). Let I = f{x) dx for some function /(x). Expressed as a first order ordinary 
differential equation, / = J(6) where dJ/dx = f{x), J{a) = 0. We solve this differential equation 
using a fourth order Runge-Kutta routine. Figure ^ demonstrates for a specific example that 
the actual numerical errors are as small as requested. All of the weak rates and thermodynamic 
quantities were calculated so that their numerical error contributions to the uncertainty in Yp were 
acceptably small. 

We compared the results of our code to the standard code, which dates back to the original 
version written in 1966 [^], was updated by Wagoner in 1973 and modernized and made 

user friendly by Kawano in 1988 |38|. Nuclear reaction rates were updated in 1993 [pl|]. One must be 
careful when making comparisons. First one must consider the numerical accuracy of the standard 
code. In 1992 Kawano [Ml estimated the accuracy of Yp to be 6 %. In 1993, Kernan addressed this 



issue in more detail and reported finding a systematic numerical error in the standard code |18, 39 1, 
6Yp = 0.0017, large enough to be very significant at our level of accuracy. Second, the standard 
code implements certain physics corrections, namely a correction put in by Wagoner to approximate 
the Coulomb correction by scaling all of the weak rates a factor, 0.98, independent of temperature. 

The systematic numerical error discovered by Kernan was measured by comparing the pre- 
dictions of the standard code at some (unspecified) integration stepsize to the predictions as the 
stepsize became very small; note, however, that the error using the default stepsize (in Ref. [^l[) is 
four times larger. The "Kernan correction" is now routinely added to the results of the standard 
code. Needless to say, a simple additive numerical correction is not adequate because other codes 
exist; not all users of the standard code use the same stepsize; and the numerical error can be 
machine dependent. 

For our comparisons we took out the Kernan and Coulomb corrections and then made the 
stepsizes small enough so that integration errors were negligible. The integration error for the 
standard stepsizes (with the two standard stepsize parameters equal to 0.3 and 0.6, respectively) 
was found to be 6Yp = 0.0073. With the standard code configured this way, we compared Yp and 
1^2/ as a function of i] in two scenarios. For the first, we used the standard weak-rate routines to 
calculate the weak rates. For the second we used our high-precision weak-rate routines to calculate 
the weak rates in the standard code. The results are shown in Fig. ^. The agreement is excellent: 
for Yp the codes differ by less than 0.15% with our weak-rate routines and by less than 0.2% with 
the standard weak-rate routines. For D the codes agree to better than 0.75%. 

This agreement gave us confidence that our code calculates Yp accurately for the baseline case 
(without the physics corrections). Of course, the convergence of two independent codes is not proof 
that they converge on the correct value. We will assume that the two codes do indeed converge 
on the the correct answer, and because our code was designed, engineered and tested for an error 
budget, we will use its results and internal error budget as the baseline for further comparison. The 
internal error budget for our code was no greater than 0.1%. 
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Figure 2: Actual numerical error in calculating Fep^un for error parameter set at ST/F = 10~^. 
The error is smaller than the specified accuracy (10~^) for all temperatures. Similar results were 
obtained for the other weak rates and thermodynamic quantities. 
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Figure 3: Comparison between the standard code and our code for ^He (lower curves) and D (upper 
curves). For the solid curves, our very accurate weak rates were inserted into the standard code. 
For the dashed curves, the standard code's weak rate routines were used. (Note: r)io = ry/lO"^^). 
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Reaction k 


6Rk/Rk 


SYp/Yp (7?io = 5.0) 


SYp/Yp (7?io = 1-8) 


n <^ p 


0.23% 


0.17% 


0.18% 




7% 


0.04% 


0.17% 


d{d, nfHe 


10% 


0.06% 


0.07% 


d{d,p)T 


10% 


0.05% 


0.06% 


Total Uncertainty 


0.19% 


0.27% 



Table 2: 1-a experimental uncertainties and their effect on Yp. All nuclear rates whose uncertain- 
ties significantly impact Yp are shown. The weak-rate uncertainty of 0.23% is due to uncertainty 
in measurements of the neutron mean lifetime, and assumes that Coulomb, radiative and ther- 
modynamic corrections to the weak rates are known to better accuracy than this. Note that for 
r] = 5.0 X 10"^'^, the neutron mean lifetime dominates the error budget. The bottom row indicates 
the RMS total uncertainty in Yp for these two values of rj. 



2.3 Nuclear Rate Uncertainties 

The primordial ^He abundance is sensitive to nuclear reactions other than the weak rates. Several 
studies of the uncertainties in theoretical abundances due to nuclear rate uncertainties have been 
performed ||l], |2l|, 39, 4C, 41, Here we will use the results and techniques of the recent work of 



Fiorentini, et al |22|. They use linear error propagation theory to quantify the effect of experimental 
uncertainties in the nuclear-reaction rates on the light element abundance uncertainties and their 
correlations, 

where the sum k is over nuclear reactions, 5Rk is the experimental uncertainty in the rate R^, and 
Afc is the logarithmic derivative 

_ log yp 

5 log Rk 



Fiorentini, et al [22| calculate the logarithmic derivatives numerically, using the standard code, and 



take the experimental rate uncertainties from Smith, et al [E^]. Contributions to the uncertainty 



in the He abundance arise almost entirely from four rates. Table 2.3 lists these rates and their 
relative experimental uncertainties. Figure ^ shows the resulting uncertainty in Yp. For rj < 
2 X 10"^'^, the reaction p{n,^)d dominates the error budget. Finally, a recent new analysis of the 
experimental uncertainties IQ, indicates the uncertainties in the reactions d{d,n)p and d{d,p)T 
have been overestimated by about a factor of two, and that the precision of the reaction p{n, 'y)d 
could be improved significantly. Thus it may well be the case that the uncertainty in Tn dominates 
the error budget for all rj. 



3 Weak Rates 

The primordial ^He abundance is very sensitive to the weak rates that maintain the balance between 
neutrons and protons. To calculate Yp to a precision of 0.12% the weak rates must be known to 
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^10 

Figure 4: The top panel shows the uncertainty in Yp due to experimental uncertainties in nuclear 
rates, as a function of r]. The solid line shows the total uncertainty, while the other lines show each 
nuclear reaction separately. The dashed line is for n p, the dashcd-dotted line is for p{n, ^)d, and 
the two dotted lines are for d{d, n)^He and d{d,p)T. The bottom panel shows the uncertainty in rj 
that would result from the above uncertainties in Yp, when rj is derived from a perfect measurement 
of the '^He abTindance. The dashed line is for the weak rate uncertainties alone, while the solid line 
is for the total nuclear rate uncertainty. The factor of ten difference in the scales between the two 
panels is indicative of the fact that Yp depends logarithmically upon rj. 
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a precision of 0.15%. In addition to numerical issues discussed earlier, several physical effects are 
important at this level: zero-temperature radiative and Coulomb corrections, finite-nucleon mass 
correction, and finite-temperature radiative correction. 

The expressions for the weak rates are derived starting with the tree-level (Born diagram) shown 
m Fig.H For purposes of illustration, we will consider the process e +p — > Ue+n. Without making 
any approximations the phase space integral for the conversion rate (per proton) can be simplified 
to a five-dimensional integral involving the matrix-element squared \A4\'^ \2G\ 

^ep^un = / dpedppd COS 9pd COS 6^d(j)y 

x|?^i \M\' fjpil - - /„), (9) 
En f^ _ (Pe + Pp) • Piy 

En V E^u 

where Eg, Ep, E^, and En denote the energies of the respective particles and J' is the Jacobian 
introduced in integrating the energy part of the delta function, and jA^I'^ is summed over initial and 
final state spins. The integration limits correspond to the kinematically allowed region in the five- 
variable phase space. An expression for E^, = py in terms of the integration variables Pe,Pp, dp, (^u, 
and (pi, is given by 



^ = i + g i- '-' -f - I. (10) 



p. 



A^B + 2E^A^ - ml{4E^ - B^) 



4^2 _ ^2 

A' = 2EeEp + - m„ — — - 2pePp cos Op, 

B = 2[p(,cos9y + Pp{cos0pcos9y + sm.9pSm9yCos4>y)\. (11) 



where E = Eg + Ep. For more details, see Ref. |26] 



This rate expression is challenging to evaluate for two reasons. First, the kinematically allowed 
region in the five-dimensional phase space is not simple. Second, the full matrix element is complex. 
Only if the nucleons are assumed to be infinitely massive, does the expression simplify: |A^|^ — > 
2^Gp{l + 3g'j^)EeEpEyEn- In that limit, the sole kinematical constraint is Ep = E^ + Q, (Q = 
nin — rup = 1.293, MeV), and the rate expression becomes a one variable integration. Normalizing 
the rates to the zero-temperature free neutron decay rate, 

i .r,„,,„(T = o)= ^?Ml±MM,„ (12) 

Ao = ^^dee(e-g)2(e2- 1)1/2 = 1.6333, ^-^3^ 

leads to the well known formula for the process ep vn: 



1 r 6(62-g2)l/2 
^^^"^^ r„Ao y„ [1 + exp(6z)] [1 + exp((g - e)z,))] ' 



(14) 



where T is the photon temperature, is the neutrino temperature, e = Ee/rue, q = Q/nie, 
z = TUe/T, and Zi, = rrie/T^- Summing the n ^ p and p ^ n rates yields the standard weak-rate 
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Figure 5: Tree level diagram for the process ep — > vn. 
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Figure 6: Weak rates as a function of temperature (Born diagram, infinite-nucleon-mass limit): (1) 
ep vn, (2) vp en, (3) en — > up, (4) vn ep, (5) n peu, (6) peu n. Note, freeze-out of 
the n/p ratio occurs at Tp ~ 0.8 MeV and ^He synthesis begins at T ~ 0.1 MeV. 



expressions H 



T„A„ V i-oo ii / (l+c«)(l + e(-i)'") ^ ' 

The six individual rates are plotted as a function of temperature in Fig. 

3.1 Zero-Temperature Coulomb and Radiative Corrections 

To order a, the weak rates with zero-temperature Coulomb and radiative corrections are given by 
the sum of the interference between the Born diagram (Fig. ^) and the diagrams in Fig. |^. 

It is conventional to separate the corrections into a Coulomb part proportional to nuclear charge 
Ze and a radiative part proportional to e. Since Z = 1 here, this separation is arbitrary. Dicus et al 



calculated the Coulomb and zero-temperature radiative corrections to the weak rates in 1982 [15|. 
Summarizing their results we obtain the following prescription for correcting the rates. First, 
perform the zero-temperature radiative corrections by multiplying the integrands of all of the rates 
by the factor, 

l + £c7(/3,y)] , (16) 



where 

-4 (2 + Up + 25/3^ + 25/3^ + 30/?^ + 20/3^ + 8p^) /(I + /3)^ (17) 
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Figure 7: Zero-temperature corrections to the process ep — ^ un. The center blob is the charged- 
current, weak- interaction vertex. 
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Figure 8: Zero-temperature radiative and Coulomb corrections to the n <-> p rates. The horizontal 
line is Wagoner's approximation to the Coulomb correction. The vertical line is at freeze-out. 



/3 is the electron's velocity and R = tanh /3 ^//3. Next apply the Coulomb correction by multiplying 
the integrand of the rates for n <-> peu and ep vn by the non-relativistic Fermi factor, 

The error from using the non-relativistic Fermi function is of order 2% of the Coulomb effect 
itself [45 1, and so the approximation is fine. Finally, Aq must be corrected for Coulomb and zero- 
temperature radiative effects by multiplying it's integrand by [l -|- ^C{P,y)] F{f3). Doing this 
increases Aq by 7.15%, to 1.7501. 

Figure ^ shows the combined zero-temperature corrections. Note that the corrections are less 
than or equal to zero for both rates for all temperatures: decreased weak rates imply earlier n/p 
freeze-out and an increase in Yp. Our code calculates the zero-temperature corrections to the 
weak rates by modifying the integrands of the rate expressions as described above, and by using 
the corrected Aq- The zero-temperature corrections yield a change, 5Yp/Yp = 1.28% which is 
insensitive to the value of rj over the range 10^1° < r] < IQ-^. This result is in agreement with 
Ref. H. 

Wagoner approximated the Coulomb correction by reducing both the n — > n and p — > n rates 
by 2%. This correction, shown by the horizontal line, is close to the high temperature asymp- 
totic Coulomb correction of —2.16%. However, n/p continues to decrease slowly for temperatures 
lower than freeze-out, where Wagoner's approximation breaks down. The fact that the real cor- 
rections are less negative in this regime means that the change in Yp from the Coulomb correc- 
tion will be less positive than one would estimate from Wagoner's approximation. Adding in the 
zero-temperature radiative corrections brings the total zero-temperature change in Yp closer to 
what would be found using Wagoner's approximation to the Coulomb correction. Table ^.l| shows 
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Correction 


5Yp/Yp 


Coulomb 
T=0 Radiative 
Combined 


1.04% 
0.24% 
1.28% 


Wagoner's approximation 


1.56% 



Table 3: Zero-temperature corrections to Yp, compared with change in Yp from Wagoner's approx- 
imation of the Coulomb correction. These corrections are insensitive to rj for lO^^^'' <f]< 10~^. 



6Yp/Yp for the Coulomb and zero-temperature radiatively separately and summed, compared to 
5Yp/Yp from Wagoner's approximation. Note in particular that the difference between Wagoner's 
approximation and the zero-temperature correction is 0.28%, which is significant at the 0.1% level. 

3.2 Finite-Nucleon Mass Correction 

Recall that the standard rate expressions, assume infinitely massive nucleons. We have 

calculated the weak rates without this assumption by numerically integrating the five-dimensional 
rate integral, Eqn.(^, using the Monte Carlo method [^]. Figure |9| shows the finite-mass corrections 
to the n <^ p rates. Using the individual rate corrections we found the corrections to the summed 
n ^ p rates. 



T Too 



n^p 



oo 



■p Too 
^ p^n J- p^n 



(19) 



(20) 



where T°° is the rate in the infinite-mass approximation, and T is the unapproximated rate. Our 
corrections are accurate to within a few percent [^6|. We incorporated the finite-mass corrections 
into our code by modifying the n <-> p rates at each temperature by the correction shown in Fig. |^. 
The resulting correction to Yp was found to be 5Yp/Yp = 0.50%, valid for 10-^° <r]< 10"^. 

3.3 Finite-Temperature Radiative Correction 

Finite-temperature modifications to the weak rates arise from several sources: 

1. the (1 lb /) quantum statistical factors in the integration over phase space 

2. a shift in the electron mass 

3. a change in the neutrino-to-photon temperature ratio 

4. a correction to the photon and fermion propagators 

5. the square of the sum of diagrams for processes that involve photons from the plasma (ab- 
sorption and stimulated emission); see Fig. |l^. 

6. finite-temperature wave-function renormalization 
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Figure 9: Finite-nucleon-mass correction to the n ^ p rates. The freeze-out temperature, Tp ~ 
0.8 MeV, is indicated with a vertical hne. 

Item H is included in our definition of the Coulomb correction. We shall define items ^ and ^ to be 
part of the thermodynamics effects, considered later. Therefore, the finite-temperature radiative 
correction to the weak rates involves items ^ |5| and ^. 

Dicus, et al [^], and Cambier, Primack and Sher |^6| calculated the finite-temperature radiative 
corrections to the weak rates. Neither of these papers correctly handle the finite-temperature wave- 
function renormalization. In fact, finite-temperature wave-function renormalization is still an open 
issue. The difficulty lies in the fact that finite temperature spoils Lorentz covariance through 
the existence of a preferred, thermal frame (in this frame the phase-space distributions are the 
Bose-Einstein or Fermi-Dirac distributions). The usual methods for obtaining the wave-function 
renormalization rely on Lorentz covariance, so that the appropriate generalization to the finite- 
temperature case is not clear. Donoghue and Holstein 17 1 start by assuming a finite-temperature 
spinor field - with creation and annihilation operators obeying the standard anti-commutation 
relations - that satisfies the nonlinear Dirac equation. They write the propagator in terms of these 
finite-temperature scalars, obtaining a finite-temperature wave-function renormalization that is a 



multiplicative factor. Sawyer |19|, and Esposito, et al [47|, start by identifying particle states with 
poles of the propagator, without reference to the finite-temperature field. They assume that the 
poles are only perturbatively shifted from their zero-temperature values. They then identify the 
finite-temperature wave-function remormalization with the residue of the propagator at the new 
pole. The result is a finite-temperature wave-function renormalization that contains additional, 
non multiplicative terms, so that the results of the two alternative approaches are different (as 
pointed out by Chapman [^]). Furthermore, the results of the Sawyer differ from Esposito, et 
al ||4^, even though they follow similar approaches. The differences change the rates for some 
processes. However, for the case of the weak rates, the three different finite-temperature wave- 
function renormalization results give the same contribution to the weak rates. For convenience, we 
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Figure 10: Finite-temperature corrections to the weak rates, i.e., corrections involving photons 
from the plasma. The bottom two diagrams represent stimulated emission. 



used the formalism of Sawyer. The correction to the en vp is given as 



5T = 
where x = m/T, 



OO fCO 



dudKpuN+{u) [N^{k^)W^{u,k.u) + Nj^{v)Wr{u,ky)] , 



(21) 



^v? - 3;2, V = v^A;2 + x"^, N±{u) = l/(e" ± 1), 



In 



U+Pu 



Wr{u,kv) 



Pu U-pu 

kyH{u) 



u-pu 
[H{u + k 



2u\ 
kv J 



[H{u + K) + H{u - K) - 2H{u)] + 



4puV 



^ Pu + ky m'^ - {uv - puky) 

2u m ; \- vln- 



Ak^PuU 



Pu ky 



and 



H(w) 



- {uv + pukv) pi - 

ij^N{-u)e{u), 
{w + q) 



(22) 
(23) 



(24) 
(25) 



with q = Q/T. The term proportional to Wr is due to finite-temperature wave function renormal- 
ization. To find the correction to the other weak rates, make the substitutions shown in Table jj. 
We calculated the finite-temperature radiative corrections to each of the weak rates. The 



correction to the summed n <-> p rates, which match Sawyer's results, are shown in Fig. 11. The 
correction formulas are complicated enough to preclude direct incorporation into our BBN code. 
Therefore we implemented these corrections as temperature-dependent fits within the BBN code. 
The resulting change in Yp, 6Yp/Yp = 0.12%, was found to be insensitive to r] in the range 
-j^g-io < ^ < 10^^. Sawyer claims a change of -|-0.02%, while Chapman claims a change of -|-0.01%. 
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Table 4: Substitutions in Eqs. (21-2^) for computing finite-temperature radiative corrections. 



Both Sawyer and Chapman compute the change in the neutron fraction to estimate 5Yp/Yp. To 
first order in the perturbation, the equations governing the evolution of the neutron fraction X„ 
and its perturbation can be written 



dT 
d5Xn 
dT 



dt 

'dT 

dt 

df 



[—Xn^n^p + (1 — Xn) F 



{5Xn + 7nXn) + Fp^„ [7^ (1 - X„) - dX^]} 



where 7„ = (5F„^p/F„^p and 7p = (5Fp^„/Fp^„. Then the change in Yp is estimated as 



6Yp SXr, 



Yp 



Xr, 



5X„ 



onset of BBN 



X„ 



(26) 



(27) 



T=0 



In order to have a direct comparison with the results of Sawyer and Chapman, we found 5Yp/Yp 



using this method. The evolution of 5Xn is shown in Fig. 12. Our results obtained from this 
approximation method confirm those using the BBN code, and differ from Sawyer and Chapman. 
However, all agree the change in Yp is small. 



4 Thermodynamics 

Thermodynamic corrections refer to corrections to the density, pressure and neutrino-to-photon 
temperature ratio. There are two effects to consider: finite-temperature QED corrections to the 
equation of state of the electromagnetic plasma, and incomplete neutrino decoupling. 

4.1 Finite-temperature QED Correction 

The finite-temperature QED corrections encompass corrections to the density, neutrino temperature 
and electron mass. All of the these corrections follow from the finite-temperature QED modification 
to the equation of state of the electromagnetic plasma. These corrections were calculated by 



Heckler |23| and applied to cosmology and solar physics. We will follow his approach, correcting a 
few small errors. 

The "^He abundance is sensitive to thermodynamic quantities in several ways. The energy density 
determines the expansion rate; changes in the expansion rate affect the freeze-out temperature, the 
abundance of free neutrons, and finally Yp. The next two effects follow from corrections to the 



18 



o 

2 4 6 8 10 



T (MeV) 

Figure 11: Finite-temperature radiative corrections to the n p rates. This plot is to be compared 
to Fig. 4 in Ref. [12]. 



X CO 



c o 



I 



.1111 I — I — I 1 1 I I I I I — I — I 1 1 I I I I I — I — I r 



'10 



I I I I I I I I I I 1_ 



0.1 



0.01 



T (MeV) 

Figure 12: Temperature evolution of the estimated change in neutron fraction Xn due to finite- 
temperature radiative corrections. The solid line shows the results of integrating the perturbation 
equations; the low-temperature asymptotic solution gives the correction to Yp, 5Yp/Yp = 5Xn/Xn- 
The arrow indicates the final result of substituting the radiative corrections into our full code. The 
two methods agree very well. 
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electron mass. A change in the electron mass affects the weak rates directly, and indirectly, by 
changing the entropy of the electron-positron plasma at the time neutrinos decouple. Since this 
entropy is transferred to the photons when the pairs disappear, this changes the neutrino-to- 
photon temperature ratio, and affects the weak rates, which are very sensitive to the neutrino 
temperature. 

The finite-temperature QED correction to the equation of state can be expressed as a modifi- 
cation to the pressure of the pressure- weighted, effective number of effective degrees of freedom, 

P{T) = Po{T)+6P{T), (28) 

where 6P(T) is the correction to the pressure and Po{T) = (vr^/OO) gpT'^ is the standard expression 
for the pressure. The change in pressure can be equated to a change in gp, dgp = 90/(7r^T^) dP. 
The correction 5P{T) can be expressed as an expansion in electron charge e ~ 0.301: 5P{T) = 
6Pi{T). The Feynman diagrams for the e^-term and e^-term are shown in Fig. 4.1. For vanishing 
chemical potential the term is 



m{T) = -— ^ / du- 



67r2 e" + 1 

' ' dudvpuPvN{u)N{v){A + In 2 , (29) 



^'^ Jx Jx \ PuPv UV - Pu Pv + 



where x = nie/T, u = Eu/T, pu = \Jv? — x"^ and N{u) = 1/(1 -|- e"). In the high-temperature limit 
T » me, 

A similar, but more involved, calculation yields the result for (5P3(T) in the limit T ^ m 



At high temperatures, the ratio 



^ ^ 11, (32) 



bP^{T) e 2 

while both the and the e^-terms are exponentially suppressed for T <^ m. Therefore, to good 
approximation, we can neglect (5^3 (T) for all T. For T ^ mf>,5gp = — 25e^/167r^. 

Prom the standard thermodynamic relation p = —P + T (dP/dT) we can find the thermody- 
namic correction to the energy density, p = po+5p, where the standard density po may be written in 
terms of the density-weighted effective number of relativistic degrees of freedom, po = (7r^/30) gp T^. 
The change in the density can be written 

30 / _ „ a _\ , 25 2 



^9P = :;^A-^P + T^SP T^^--^^e\ (33) 



Figure 14 shows 6gp and 6gp as a function of temperature. 
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Figure 14: Finite-temperature QED change in pressure-weighted {gp, soUd hne) and density- 
weighted {gp, dashed hne) relativistic degrees of freedom. 
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Figure 15: The top panel shows the finite-temperature QED correction to the electron mass as a 
function of temperature. The dashed curve neglects the p-dependent term, while the sohd curve 
assumes p = 3T. The bottom panel shows the relative error due to not including the p-dependent 
term. This error, which is a ten percent correction to the correction, can be safely neglected. 



The finite-temperature QED correction to the pressure is a change in the dispersion relation of 
the electrons which can be attributed to a change in the electron mass: 



■ + w? + 5vn? . 



(34) 



The formula for 5m? follows from the definition of the pressure correction |48|. 

1 



6m {p, T) 



6 

27r2p 



oo r, 
, run 

du 

u + 1 



du In 



Pu + 



Pu 



1 



(35) 



and Pu = p/T. Figure 15 shows the finite-temperature QED 



where X — m^^jT^ ku 

correction to the electron mass as a function of temperature. Figure 16 shows the effect of the shift 
in the electron's mass on the n <-> p rates. The lower curves indicate the error due to not including 
the momentum-dependent part of the mass correction. For our calculations, the error is negligible 
and we neglect the p-dependent term in the mass correction formula. 

The final effect of the thermodynamic corrections is a change in the neutrino-to-photon temper- 
ature ratio. This can be derived starting with the expression for 5P{T) and tracking the entropy 
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T (MeV) 

Figure 16: The top curves show the effect of the finite-temperature electron-mass correction on the 
weak rates. The sohd curve is for n ^ p and the dashed curve is for p ^ n. The bottom curves 
show the error due to not including the p-dependent term in the mass correction formula. 



density of the neutrinos and other particles. Let Si, be the entropy density of neutrinos and sem 
be the combined entropy density of the electrons, positrons and photons: 



•Sem 



Pu + Pu ^ TtT^ 3 

30 " ' 

Pe± + Pe± +P^ + 



47r2 _2_ 
15" ^ 37r2 



oo / 2 _ J.2 

du ^ (4^ 

e" + 1 ^ 



vr 



u 



) + ^ (Sgp + ^Sgp 



(36) 



(37) 



In the limit that the neutrinos are completely decoupled, the two entropies per comoving volume 
are separately conserved: SuO,^, sem^^ = constant, where a is the scale factor. The small residual 
coupling of the neutrinos to the electromagnetic plasma leads to a correction of about ~ 0.1% [27|, 
discussed below, which can be ignored here. At high temperature we have 



sem« 



while for all temperatures, 



| + ^l*»p(r) + 3%,(r)l.| 



1 



25 e2 
88^ 



(38) 



n 



21 ^ 7^ 



du 



e" -Fl ^ 



)+-[5gp{T) + 35gp{T)] 



(39) 



Assuming that the neutrinos decouple at a temperature To ~ 2 MeV ^ mg and taking the ratio of 
entropies to be given by Eqn. (Ssl), it follows that the ratio of the neutrino-to-photon temperature 
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Figure 17: Relative finite-temperature QED change in the neutrino temperature, as a function of 
photon temperature. Note that the zero-temperature Umit is altered from the standard value by 
about 0.08 %. 
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The zero-temperature limit of the neutrino temperature photon temperature relation is altered 
This makes sense physically: the positive correction to the electron mass means that the electron- 
positron plasma has less entropy to give to the photons upon annihilation, and thus photons are 
heated less than they would be without the correction. Figure |l^ shows the finite-temperature 
QED change in neutrino temperature versus photon temperature. 

We incorporated the QED corrections to the equation of state into our code by changing the 
energy density, the electron mass in the weak-rate calculations and the neutrino temperature. 
The resulting change in Yp, 5Yp/Yp = +0.043% was found to be insensitive to r] in the range, 
10^^'^ < 7] < 10~^. Dicus, et al |15| attempted to calculate the thermodynamic corrections, and 
found 6Yp/Yp = —0.04%, but only included the effect of the electron mass on the weak rates. 
Heckler estimated the effect on Yp and found 5Yp/Yp = +0.06%. (It should be noted that his 
value for the change in neutrino temperature was incorrect.) In any event, the thermodynamic 
correction to Yp is small. 



This expression differs somewhat from the result obtained by Heckler |23(]. He now agrees with our result 
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4.2 Incomplete Neutrino Decoupling 

The standard code assumes that neutrinos decoupled completely before annihilations. It has 
been pointed out that this assumption is not strictly valid ||T^. Neutrinos are "slightly coupled" 
when pairs are annihilated, and hence share somewhat in the heat released. The first calculations 
15, E^, 5C] of this effect were "one-zone" estimates that evolved integrated quantities through the 



process of neutrino decoupling. More refined "multi-zone" calculations tracked many energy bins. 



assumed Boltzmann statistics and made other approximations. |27, 51 1. The latest refinements 



have included these small effects as well [q^, |53|, gj]. Fields et al ||55|] incorporated the slight effect 



of the heating of neutrinos by annihilations into the standard code and found a shift in ^He 
production, 5Yp = +1.5 x 10"'*, which is insensitive to rj for 10"^'' l^f] < 10~^. 

5 Summary 

All of the physics corrections we investigated have been studied elsewhere. However, not all of 
them have been implemented in a full code; some have been implemented incorrectly; and there 
have been changes in some of the physics corrections. Further, the issue of numerical accuracy of 
the standard code has not been comprehensively and coherently addressed. Finally, the corrections 
have been implemented in a patchwork fashion, so that the users of many codes do not know which 
corrections are in, which are out, and which may be double counted (e.g., by adding the numerical 
correction and running a small stepsize). As noted earlier, the results of a number of BBN codes 
gave a 1% spread in the prediction for Yp with the same value of r] and r„. 

The goal of this work was a calculation of the primordial ^He abundance to a precision lim- 
ited by the uncertainty in the neutron mean lifetime, (5r„ = ±2 sec, or 5Yp/Yp ~ 0.2%, with 
reliable estimates of the theoretical error. To achieve this goal we created a new BBN code, de- 
signed, engineered and tested to this numerical accuracy. To this baseline code we added the 
microphysics necessary to achieve our accuracy goal - Coulomb and zero-temperature radiative 
corrections, finite-nucleon-mass corrections, finite-temperature radiative corrections, QED thermo- 
dynamical corrections, and the slight heating of neutrinos by annihilations. These corrections - 
coincidentally all positive - increase the predicted ^He abundance by 5Yp = 0.0049 or 2%. Table ^ 
summarizes these corrections for rj = 5x 10^^*^. For each physical or numerical effect, we have been 
careful to control the error in Yp introduced by approximations or inaccuracies to be well below 
0.1%. With confidence we can state that the total theoretical uncertainty is less than 0.1%. 

Summarizing our work in one number 

Yp{r] = 5 X 10"^°) = 0.2462 ± 0.0004 (expt) ± < 0.0002 (theory). (42) 
Further, the precise value of the baryon density inferred from the Buries- Tytler determination 



of the primordial D abundance, O^^^ = 0.019 it 0.001 |^, 57 1, leads to the prediction: Yp = 
0.2464 ± 0.0004 (expt) ± 0.0005 (D/H) ± < 0.0002 (theory). 

Finally, we give two fitting formulae for our high-accuracy ^He predictions. The first, is accurate 
to better than 0.05% and is valid for 10"^° < r] < 10"^, = 3.00 and 880 sec < r„ < 890 sec. In 
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Yp 


Cumulative 
5Yp(x 10-4) dYp/Yp{%) 


Effect Alone 
(5yp(x 10-4) SYp/Yp{%) 


Baseline 


0.2414 










Coulomb and T = radiative 


0.2445 


+31 


+1.28 


+31 


+1.28 


finite mass 


0.2457 


+43 


+1.78 


+12 


+0.50 


finite T radiative 


0.2460 


+46 


+1.90 


+3 


+0.12 


QED plasma 


0.2461 


+47 


+1.94 


+1 


+0.04 


residual z^-heating 


0.2462 


+49 


+2.00 


+1.5 


+0.06 



Table 5: Summary of results. For absolute numbers we have picked t] = 5.0 x 10- . By baseline 
we mean the results of our BBN code without any of the physics effects listed, and with small 
numerical errors. 



terms of C = 10 + log^o ??, 



yp(C,r„) = yp (C, 885.4 sec) + (rn- 885.4 sec) (5Yp(C), 
Yp(C, 885.4 sec) = (ao + ai C + a2 + ag + 04 C^) , 

SYpiO = (6o + 6iC + ?'2C' + &3C' + &4C') (43) 



where the coefficients Oj, bi are given by 



ao = 0.22292 , oi = 0.05547 , as = -0.05639 , 

03 = 0.04587 , 04 = -0.001501 

60 = 2.082 X 10-4 , 61 = 0.535 x 10-^ , 62 = -2.856 x 10-^ , 

63 = 4.672 X 10-4 , 64 = 2.420 x 10-4 . 



(44) 



The second fitting formula is accurate to 0.5% and is valid for 10 < 77 < 10 ^, 880 sec < < 
890 sec, and 2.5 < A^^ < 4.0. 

Yp{C, T, N,) = Yp{C, T, 3) + {N, - 3) (co + ci C + C2 C' + C3 + C4 C^) , (45) 

where 

Co = 0.01276 , ci = 0.00409 , C2 = -0.00703 , 

C3 = 0.00571 , C4 = -0.00186 . ^ ' 

Acknowledgements 

The BBN code we described in this paper began as a project for a graduate course in cosmology. 
R.L. gratefully acknowledges the work of his Red Team colleagues, Moses Hohman, Mike Joffre, 
Russel Strickland and Craig Wiegert. Thanks to Gary Steigman for helpful conversations. This 
work was supported in part by the DOE (at Chicago and Fermilab) and by NASA at Fermilab 
through grant NAG 5-2788, and 5-7092. 



26 



References 

[1] D.N. Schramm and M.S. Turner, Rev. Mod. Phys. 70, 303 (1998). 

[2] J. Yang, M.S. Turner, G. Steigman, D.N. Schramm, and K.A. Ohve, Astrophys. J. 281, 493 
(1984). 

[3] T.P. Walker et al., Astrophys. J. 376, 51 (1991). 

[4] Y. Izotov and T.X. Thuan, Astrophys. J. 500, 188 (1998). 

[5] Y. Izotov and T.X. Thuan, Astrophys. J. 497, 227 (1998). 

[6] K. Ohve and G. Steigman, Astrophys. J. Suppl. 97, 49 (1995). 

[7] K. Olive, E. Skihman, and G. Steigman, Astrophys. J. Suppl. 483, 788 (1997). 

[8] K. Ohve, E. SkiUman, and G. Steigman, Astrophys. J. 489, 1006 (1997). 

[9] B.E.J. Pagel and A. Kazlauskas, MNRAS 256, 49 (1992). 
[10] B.E.J. Pagel et al., MNRAS 255, 325 (1992). 
[11] S. Sarkar, Rep. Prog. Phys. 59, 1493 (1996). 

[12] C.J. Copi, M.S. Turner, and D.N. Schramm, Phys. Rev. D 55 (1997). 
[13] V.F. Schvartsman, JETP Lett. 9, 154 (1969). 

[14] G. Steigman, D.N. Schramm, and J. Gunn, Phys. Lett. 66B, 202 (1977). 

[15] D. Dicus et al., Phys. Rev. D 26, 2694 (1982). 

[16] J.F. Donoghue and B.R. Holstein, Phys. Rev. D 28, 340 (1983). 

[17] J.F. Donoghue and B.R. Holstein, Phys. Rev. D 29, 3004 (1984). 

[18] P. Kernan, PhD thesis, Ohio State University, 1993. 

[19] R.F. Sawyer, Phys. Rev. D 53, 4232 (1996). 

[20] LA. Chapman, Phys. Rev. D 55, 6287 (1997). 

[21] M.S. Smith, L.H. Kawano, and R.A. Malaney, Astrophys. J. Suppl. 85, 219 (1993). 

[22] G. Fiorentini, E. Lisi, S. Sarkar, and F. L. Villante, Phys. Rev. D 58, 063506 (1998). 

[23] A.J. Heckler, Phys. Rev. D 49, 611 (1994). 

[24] D. Seckel, hep-ph/9305311 (1993). 

[25] G. Gyuk and M.S. Turner, astro-ph/9307025 (1993). 

27 



[26] R.E. Lopez, M.S. Turner, and G. Gyuk, Phys. Rev. D 56 (1997). 

[27] S. Dodclson and M.S. Turner, Phys. Rev. D 46, 3372 (1992). 
[28] C. Caso et al., Euro. Phys. J. C3, 1 (1998). 

[29] S. Arzumanov et al., in V International Seminar on Interaction of Neutrons with Nuclei, 
Dubna, 1997. 

[30] J. Byrne et al., Europhys. Lett. 33, 187 (1996). 

[31] L. Kawano, FermiIab-pub-92/04-a, 1992. 

[32] N. Itoh, A. Nishikawa, S. Nozawa, and Y. Kohyama, Astrophys. J. 488, 507 (1997). 
[33] R.V. Wagoner, W. Fowler, and F. Hoyle, Astrophys. J. 148, 3 (1966). 
[34] R.V. Wagoner, Astrophys. J. 179, 343 (1973). 
[35] J. Lara, astro-ph/9806057 . 

[36] W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, Numerical Recipes in C, 
Cambridge University Press, Cambridge, 1990. 

[37] R.V. Wagoner and D.N. Schramm, Ann. Rev. Nucl. Part. Sci. 27, 37 (1977). 

[38] L. Kawano, Fermilab-pub-88/34-a, 1988. 

[39] P.J. Kernan and L.M. Krauss, Phys. Rev. Lett. 72, 3309 (1994). 

[40] L. Krauss and P. Romanelh, Astrophys. J. 358, 47 (1990). 

[41] L.M. Krauss and P.J. Kernan, Astrophys. J. Lett. 432, 79 (1994). 

[42] N. Hata et al., Phys. Rev. Lett. 75, 3977 (1995). 

[43] S. Buries, K. Nollett, J. Truran, and M.S. Turner, astro-ph/9901157 (1999). 

[44] S. Weinberg, Gravitation and Cosmology, J. Wiley, New York, 1972. 

[45] D.H. Wilkinson, Nucl. Phys. A 377, 474 (1982). 

[46] J. Cambier, J.R. Primack, and M. Sher, Nucl. Phys. B 209, 372 (1982). 

[47] S. Esposito, G. Mangano, G. Miele, and O. Pisanti, 1998. 

[48] J.I. Kapusta, Finite- Temperature Quantum Field Theory, Cambridge University Press, Cam- 
bridge, 1989. 

[49] M. A. Herrera and S. Hacyan, Astrophys. J. 336, 539 (1989). 

[50] N. C. Rana and B. Mitra, Phys. Rev. D 44, 393 (1991). 

28 



[51] A. D. Dolgov and M. Fukugita, Phys. Rev. D 46, 5378 (1992). 

[52] N.Y. Gnedin and O.Y. Gnedin, Astrophys. J. 509, 11 (1998). 

[53] S. Hannestad and J. Madsen, Phys. Rev. D 52, 1764 (1995). 

[54] A. D. Dolgov, S. H. Hansen, and D.V. Semikoz, Nucl. Phys. B 503, 426 (1997). 

[55] B. Fields, S. Dodelson, and M.S. Turner, Phys. Rev. D 47, 4309 (1993). 

[56] S. Buries and D. Tytler, Astrophys. J. 499, 699 (1998). 

[57] S. Buries and D. Tytler, Astrophys. J. 507, 732 (1998). 



29 



